# C1
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

slope_temp <- 'separate'
order_temp <- 1

dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg1 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)

contVar <- cont1
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg2 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont2
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg3 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont3
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont4
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg5 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont5
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg6 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont6
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg7 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont7
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg8 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

fileName <- paste0('results/APPENDIX_C/C1.html')
reg <- list(reg1, reg2, reg3, reg4, reg5, reg6, reg7, reg8)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# C2: UPDATE!
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

contVar <- contVarBase
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))

reg1 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 1, covariates = contVar %>% paste(., collapse = '+'))
reg2 <- dat_temp1 %>% rdd_reg_np(slope = 'separate', bw = bw, covariates = contVar %>% paste(., collapse = '+'))
reg3 <- dat_temp1 %>% rdd_reg_np(slope = 'same', bw = bw, covariates = contVar %>% paste(., collapse = '+'))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 2, covariates = contVar %>% paste(., collapse = '+'))
reg5 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 2, covariates = contVar %>% paste(., collapse = '+'))
reg6 <- dat_temp1 %>% rdd_gen_reg(slope = 'separate', bw = bw, order = 1, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))
reg7 <- dat_temp1 %>% rdd_gen_reg(slope = 'same', bw = bw, order = 2, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))
reg8 <- dat_temp1 %>% rdd_gen_reg(slope = 'separate', bw = bw, order = 2, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))

# LINEAR
fileName <- paste0('results/APPENDIX_C/C2_1.html')
reg <- list(reg1, reg4, reg5)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

# NON-PARAMETRIC CANNOT BE EXPORTED

# LOGIT
fileName <- paste0('results/APPENDIX_C/C2_2.html')
reg <- list(reg6, reg7, reg8)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = T,
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# C3: UPDATE
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

dat_temp1 <- rdd_data(df$y, df$x, covar = df[,c('pid', contVar)], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg1 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 1, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg2 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 1, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg3 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 2, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 2, 
                                 covariates = contVar %>% paste(., collapse = '+'))

fileName <- paste0('results/APPENDIX_C/C3.html')
reg <- list(reg1, reg2, reg3, reg4) %>% lapply(., function(x) as.lm)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
reg <- list(reg1, reg2, reg3, reg4)
modelsummary(reg, 
             vcov = list(vcovCL_custom(reg1), vcovCL_custom(reg2), vcovCL_custom(reg3), vcovCL_custom(reg4)), 
             exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# C4-C6
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

slope_temp <- "separate"
order_temp <- 1

# LINEAR
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = month_seq[b], order = order_temp, 
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]]$param <- month_seq[b]
}

fileName <- paste0('results/APPENDIX_C/C5&6.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  createATTPlotBW(., ybreaks = seq(-2.5, 2.5, 0.05), ylimits = c(-0.1, 0.1), yintercept = 0)
fileName_plot <- paste0('results/APPENDIX_C/C4_A.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

# NONPARAMETRIC
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_reg_np(slope = slope_temp, bw = month_seq[b], 
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]] <- reg[[b]][['RDDslot']][['model']]
  reg[[b]]$param <- month_seq[b]
}

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  createATTPlotBW(., ybreaks = seq(-2.5, 2.5, 0.05), ylimits = c(-0.1, 0.1), yintercept = 0)
fileName_plot <- paste0('results/APPENDIX_C/C4_B.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

# LOGIT
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = month_seq[b], order = order_temp, 
                                        fun = glm, family = binomial(link = 'logit'),
                                        covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]]$param <- month_seq[b]
}

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = T) %>%
  createATTPlotBW(., ybreaks = seq(-5.5, 5.5, 0.25), ylimits = c(0.5, 2), yintercept = 1)
fileName_plot <- paste0('results/APPENDIX_C/C4_C.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

############

# C7
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

mediatorvec <- c("age_olderthan47", "nchild_lessthan2", "nchild_morethan1", 
                 "married_married", "married_notmarried", "income2_lessthan850")
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_C/C7.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# C8
source('scripts/createDataset.R')
df <- df_turnout %>% mutate(y = turnout13)

mediatorvec <- c("age_olderthan47", "nchild_lessthan2", "nchild_morethan1", 
                 "married_married", "married_notmarried", "income2_lessthan850")
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = bw, order = order_temp,
                                          fun = glm, family = binomial(link = 'logit'))
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = bw, order = order_temp,
                                        fun = glm, family = binomial(link = 'logit'))
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_C/C8.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = T, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)
